% This script finds the following for every Rz:
% - perceived utility from a contract given f = 0;
% - a maximum fee fmax that an agent is willing to pay for a given return,
% - a maximum EARLY fee that is calculated under the assumption that fees
% are collected only until retirement. 

clear
load('recalibrated_data.mat')

% 1) This part calculates the lifetime utility for any contract (as if f=0)
rho = mydata.rho;
delta = mydata.delta;
effhh_size_ = mydata.effhh_size;
avgIncome = mydata.avgIncome;
meanhhs = mydata.meanhhs;
R = mydata.R;
zpen_eval = mydata.zpen_eval;
alpha = mydata.alpha;
survival = mydata.survival;

% 1a) The perceived utility from the outside option

%instantaneous utility from consumption at each point in time
instant = zeros(1,71);
for t = 1:71
    if rho == 1
        instant(t) = effhh_size_(t) * log(mydata.outside_avgTC_TC(t)/effhh_size_(t));
    else
        instant(t) = effhh_size_(t) * ((mydata.outside_avgTC_TC(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
    end
end

% utility for leaving a bequest at each point in time 
if rho == 1
    yearly_nobeq = meanhhs * log(avgIncome/meanhhs);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
else
    yearly_nobeq = meanhhs * ((avgIncome / meanhhs)^(1-rho)-1) / (1-rho);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
end

bequest = zeros(1,71);
for t = 1:71
    cons_bequest = (R-1)*(max(0,mydata.outside_avgX_TC(t)+(1-zpen_eval(t))*mydata.outside_avgZ_TC(t))) + avgIncome;
    if rho == 1
        yearly_withbeq = meanhhs * log(cons_bequest/meanhhs);
        utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
    else
        yearly_withbeq = meanhhs * ((cons_bequest / meanhhs)^(1-rho)-1) / (1-rho);
        utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
    end
    bequest(t) = alpha*(utility_withbeq - utility_nobeq);
end

% put together lifetime utility
lifetime = instant(1);
for t = 2:71
    lifetime = lifetime + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
end

outside_lifetime = lifetime;
mydata.outside_lifetime = outside_lifetime;

% 1b) The perceived utility for any Rz and fee = 0, the maximum fee willing
% to pay for any Rz, the maximim early fee willing to pay.
for THETA = mydata.THETA_grid
    suffix = strcat('avgZ',num2str(THETA*10000),'_TC');
    avgZ_TC = mydata.(suffix);
    suffix = strcat('avgX',num2str(THETA*10000),'_TC');
    avgX_TC = mydata.(suffix);
    suffix = strcat('avgTC',num2str(THETA*10000),'_TC');
    avgTC_TC = mydata.(suffix);
    
    instant = zeros(1,71);
    for t = 1:71
        if rho == 1
            instant(t) = effhh_size_(t) * log(avgTC_TC(t)/effhh_size_(t));
        else
            instant(t) = effhh_size_(t) * ((avgTC_TC(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
        end
    end
    
    bequest = zeros(1,71);
    for t = 1:71
        cons_bequest = (R-1)*(max(0,avgX_TC(t)+(1-zpen_eval(t))*avgZ_TC(t))) + avgIncome;
        if rho == 1
            yearly_withbeq = meanhhs * log(cons_bequest/meanhhs);
            utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
        else
            yearly_withbeq = meanhhs * ((cons_bequest / meanhhs)^(1-rho)-1) / (1-rho);
            utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
        end
        bequest(t) = alpha*(utility_withbeq - utility_nobeq);
    end
    
    lifetime = instant(1);
    for t = 2:71
        lifetime = lifetime + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
    end
    
    suffix = strcat('lifetime_',num2str(THETA*10000));
    mydata.(suffix) = lifetime;
    
    % Now onto the maximum fee
    fmax = 0;
    lifetime_fee = lifetime;
    while lifetime_fee > outside_lifetime
        avgTC_fee = avgTC_TC - fmax*ones(1,71);
        
        instant = zeros(1,71);
        for t = 1:71
            if rho == 1
            instant(t) = effhh_size_(t) * log(avgTC_fee(t)/effhh_size_(t));
            else
            instant(t) = effhh_size_(t) * ((avgTC_fee(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
            end
        end
    
        lifetime_fee = instant(1);
        for t = 2:71
            lifetime_fee = lifetime_fee + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
        end
        
        fmax = fmax +10;
    end
    
    suffix = strcat('fmax_',num2str(THETA*10000));
    mydata.(suffix) = fmax - 15;
    
    % And finally, maximum fee when fees are charged PRIOR TO RETIREMENT
    % only
    fmax_early = 0;
    lifetime_fee_early = lifetime;
    while lifetime_fee_early > outside_lifetime
        avgTC_fee_early = avgTC_TC - fmax_early*[ones(1,45),zeros(1,71-45)];
        
        instant = zeros(1,71);
        for t = 1:71
            if rho == 1
            instant(t) = effhh_size_(t) * log(avgTC_fee_early(t)/effhh_size_(t));
            else
            instant(t) = effhh_size_(t) * ((avgTC_fee_early(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
            end
        end
    
        lifetime_fee_early = instant(1);
        for t = 2:71
            lifetime_fee_early = lifetime_fee_early + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
        end
        
        fmax_early = fmax_early +10;
    end
    
    suffix = strcat('fmax_early_',num2str(THETA*10000));
    mydata.(suffix) = fmax_early - 15;
end


save('recalibrated_data_with_eval','mydata')

    
    

